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ABSTRACT 


When the signal to noise ratio is low, the detection of a track in a lofargram by computer 
is dinicult. Actually, operators often face this kind of problem when detecting an 
acoustic signature in the noisy sea environment. Furthermore, operators must keep on 
tracking lofargram to identify a given target. This problem could be handled by the au¬ 
tomation of lofargram processing using filtering or image processing techniques. This 
technique will suppress the background and emphasize the spectral lines in the 
lofargram. Targets can be tracked by an automatic lofargram processing system up to a 
certain point at which the system should alarm the operator. The enliancement proc¬ 
essing of lofargrams using a relaxation method could be one part of the automatic sys¬ 
tem of lofargram processing to provide available target information for good decisions. 
The objective of this thesis is to enhance spectral lines of the lofargram by using the re¬ 
laxation method which is an iterative approach to line detection. This technique makes 
probabilistic decisions at every point of the lofargram for each iteration. Decisions are 
adjusted at sur‘ ''ssive iterations based on the results of the previous iterations. This 
thesis tested algorithms using relaxation method for lofargrams. Some experimental re¬ 
sults of the lofargram processing are presented. Experimental results showed that an 
edge relaxation method can '.'ield better results than the line relaxation method. How¬ 
ever, a double line detected from an edge is still undesirable and requires future work. 
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1. INTRODUCTION 


A. WHAT IS A LOFAR 

The term "LOFAR" is flefined as "search technique using omnidirectional 
sonobuoys". It is an acronymn for "LOw Frequency Analysis and Recording." Every 
ship including submarine create peculiar noises in the sea. The noise spectrum ran be 
used to distinguish dilTcrent types of snips. If we know the noise frequency spectrum 
pattern radiated from a certain ship, it i^ possible to identify this ship by examining the 
lofargram. 'I'his is the basic motive of using lofar. 

For determining the characteristic of noise spectra, a frequency-time analyzer can 
be used, which is similar to those used for speech analysis. It yields a plot of frequency 
against time and shows the intensity of the sound in the bandwidth of the spectrum by 
darkening the record paper. This paper plot is called a lofargram. 

Before the appearance of the digital lofargram, thermal burning of record paper was 
used. But, the conventional paper gram was inconvenient to handle compared to digital 
lofargram. Furthermore, an undesirable smell was created during operation because of 
the burning on the long strip of paper using the old method. Therefore, the old method 
is almost replaced by the digital lofargram. The digital lofargram is obtained by digitizing 
the analog signal and displaying a gray .scale on the computer screen or on a paper. 

Figure 1 adopted from [Ref. 1] is an example of a lofargram of the noise of a large 
surface ship obtained from a hydrophone. As the ship passed by, tracks are left in the 
lofargram. The frequency scale ,dong the horizontal axis extends from 0 to 150 Hz anci 
the recording duration of the vertical axis was approximately 1/2 hour. The constantly 
spaced vertical line components marked by arrows aie blade-rate lines. The linos marked 
X are of unknown origin. Figure 2 adopted from [Ref. 2] indicates the block diagram of 
a passive sonar to get the lofargram. 

Generally, the lofargram has characteristics as follows 

1. Signal to noise ratio is generally low. 

2. A constant frequency tonal produces a darkening over time periods, which would 
appear as a vertical line on the lofargiam display. 

3. Multi-lines could be displayed on the lofargram depending on the spectrum band¬ 
width. 

4. Slant lines or slightly curved lines arc due to the Doppler shift could be shown on 
the lofargram. 
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Figure 1. Lofargrani of a surface ship at a speed of 11 knots 



Figure 2. Passive sonar system block diagram 


B. OBJECTIVE OF LOFARGRAM PROCESSING 

When the SNR is lower than 3db. the interpretation of lofargrani is difficult. Actu¬ 
ally, operators often face this kind of situation in the noisy sea environment. Further¬ 
more, operators must keep on tracking lofargram to identify target. This problem could 
be handled by the automation of lofargram processing using filtering or image process¬ 
ing techniques which can suppicss the background and emphasize the spectral lines in 
the lofargram. Target can be tracked by an automatic lofargram processing system up 















to a certain point then the system alarms the operator. The enhancement processing of 
a lofargram using the rela.xation method could be one part of the automatic system of 
lofargram processing to provide target related information for good decisions. 

The objective of this thesis is to enhance the spectral lines of the lofargram by using 
the relaxation method which is an iterative approach to line detection. This technique 
makes probabilistic decisions at every point of the lofargram for each iteration. Deci¬ 
sions are adjusted at successive iterations based on results of the previous iterations. 

This thesis will present some experimental results of the line enhanced lofargrams 
using the relaxation method. 

C. WHAT IS EDGE AND LINE DETECTION 

Edges are defined as the discontinuity of intensity in an image, and arc basic parts 
of the image information in general. Both the edge and line detection are fundamental 
techniques of image processing which are used to detect boundary of an object. But, 
there is an important difference between edge and line detection. Edge detection aims 
at finding the local discontinuities in an image. These discontinuities arc of interests be¬ 
cause they are likely to occur at the boundaries of objects. The output of an edge de¬ 
tection process will be some local line shown as boundaries of an object, figure 3 shows 
the step edge and output profile of an edge detection process. 



Figure 3. Gray lesel profile of a) edge and b) line 

In real images, with noise and surface imperfections, gaps between local lines arc ex¬ 
pected. Therefore, a process to organize the local edges into aggregates is needed in a 
line detection technique. Line detection is a process where edges in an image are ag¬ 
gregated to form object boundaries. These boundaries may not be known a priori. But, 
in many cases, they can be appro,ximated well by piecewise linear .segments. Mouever, 
it is not feasible to simply fit linear segments to all the edges in an image and discard the 
poor fits. It is first necessary to aggregate the edges lying along a single line or along 
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another known curve. Proximity and the directions of an edge or more detailed de¬ 
scriptions of the edges can be used for such aggregation. Therefore, edge detection and 
line detection can be put in sequence. Basic ideas of edge and line detection are ex¬ 
plained in (Ref. 3]. 

Usually edge and line detection algorithms can be divided into the following cate¬ 
gories. [Ref. 4] 

1. Detection of element. 

Difierential-type operator. 

Model-fitting method. 

2. Enhancement of lines. 

Etthancement in image space. 

Enhancement in parametric space (Hough transform). 

3. Connection of elements. 

Connections using graph search techniques. 

Category 1 refers to the process to detect line or edge elements by applying local oper¬ 
ators such as Roberts (Ref. 5j , Prewitt (Ref. 6] and Sobel (Ref. 7J. An approximation 
of the gradient for a digital picture from Roberts's opciator is given by 

R{iJ) =•• 

- V{«('+ i)-s(v)}V(g(/,./+ I)-g(/+ lj)f (1.1) 


In Figure 4 g{ij) denotee the intensity of the itnage iof the pixel (i,j), which the direction 
of the gradient is given by 


<? = - 


_7t 

4 


+ tair’[ 


g{ij+ l)-g(<’+ li/) 
g{i +IJ+ i)-g(v) 


] 


( 1 . 2 ) 


'fhe Prewitt operator defines the magnitude of a gradient from Figure 5 by 

5-7^+^ 


(1.3) 


Sx = (^ 2+034 a^yia^+a^+Of,) 


(1.4) 


Sy = {a(,■\-aya^)~{aya^+a^) 


(1.5) 


and the direction of the edge is given by 
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(1.6) 
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« 


6 = tan '( 



) 



Figure 4. Robert's gradient operator 

The Sobel operator delines a gradient from figure 5 by 

S = I (c:y+2u,+<72)-(Ot,+2as+(74) | 

+ I (uo+2a/4-a6)-(rt2+2o3-^-c74) I (1.7) ' 

Generally, a differential type operator doesn't perform well compared to a model fitting 
method. We are going to discuss this problem in the e.xperiment part of this thesis. 

Category 2 is a smoothing processes which eliminates noise components and em¬ 
phasize edges. In this category, the iterative method performs edge enhancement by ap¬ 
plying iterative processing based on relationships of its neighbour elements in an image 
space. Category 3 is a process which connects small line segments into longer lines. 
Knowledge or constraints related to the sequence of elements to be detected are used in 
an evaluation function. Then, the connection of elements is regarded as a graph search 
problem. In this case, two techniques are used: one directly proceeds with the line ex¬ 
traction while searching the edge elements for brightness data, and the other applies a 
local edge operator to the entire image before the connection of elements is done. 
[Ref.3] 
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Figure 5. Apixel and Its eight neighbors 

In this chapter, the general ideas of line and edge detection as used in lofargratn 
processing arc explained. Chapter II introduces the relaxation technique, while a modi¬ 
fied algorithms using relaxation technique is introduced in Chapter III. Experimental 
results of lofargram processing are given in Chapter IV. Chapter V is the conclusion 

of this thesis. 
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II. RELAXATION METHOD 


A. GENERAL BACKGROUND 

Relaxation tccliniqucs were first used by Waltz [Ref. 8] for the description of solids 
and were expanded by Zuckcr [Ref. 9] for a variety of applications. Relaxation is an it¬ 
erative approach to assign each pixel to categories by assigning the neighboring pixels 
in a "compatible" way. 1 his means that each pixel has interaction with its neighboring 
pixels, and the degree of interaction is considered in terms of the compatibility of each 
pixel to its neighboring pixels, fhis method takes a set of probabilities for each pixel to 
belong to a possible class, and uses an iteiative technique to update the probabilities. 

An important clement of the method is a set of compatibility measures c(/<,.s;/,/) that 
gives the compatibility of assigning pixel k to class s and assigning the neighboring pixel 
/ to class /. It is assumed that c(/<,s;/,/) is in the range (-1,+ 1) where -1 implies a strong 
incompatibility, + 1 implies a strong compatibility, and 0 implies neutrality. p'^{k,s) is an 
initial estimate of the probability that pixel k belongs to class s. Let q'{k,3) is the incre¬ 
ment of the probability and N bt the itumber of neighbors for each pixel; a neighbor¬ 
hood of 3x3 or 5x5 is usually employed. Then, a sequence of estimates p'{k,s) is 
computed iteratively as follows 

1. For each pixel /<, compute the "neighbor compatibility" between pixel k t.nd its 
neighbor pixel / in class s, which results in the probability of pixel k in class s. 

= ( 2 . 1 ) 

/ f 

where the summation is over all possible neighbor pixels / and all classes t. If only 
two classes(for example, edge and non-edge) are considered, then there are only 
four cases of neighbor compatibility c{k,e\l,e) means the compatibility of an edge 
pixel k being a neighbor of an edge pixel /. c{k,c:,l,n) means the compatibility of an 
edge pixel k being a neighbor of a non-edge pixel I . c{k,n;l,e) means the compat¬ 
ibility of a non-edge pixel k being a neighbor of an edge pixel /. c{k,iv,l,n) means 
the compatibility of a non-edge pixel /< being a neighbor of a non-edge pixel / . 
Therfore, equation (2.1) can be expanded as follows 
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( 2 . 2 ) 


q\k,e) = — ^c(/c, <?;/,«)//(/, e)+c(/<,e;/,H)/?V.«) 
/ 


q'{k,n) = ~ '^c{k,n;l,c)p^(l,e)+c{k,n\l,n)p{l,it) (2.3) 

i 


2. Update the probabilities p'*'{k,s) 

r+\,, _ r\k,s)iiWii<,s)) 

P \ktS) 

2_^p\k,s)[l+q\k,s)) 

S 


Vor two classes; 

r+1.^ . ^_ />V<,g)[l+?Vc,c)] _ 

^ />V<,e)Li+'/'(/<.c)J+/>'(^.«)[J+<7'(^.»)] 


(2.4) 


(2.5) 


Equauon (2.5) gives the probability that pixel k is an edge. Then the probability 
that pixel /< is a non-edge is defined by 

//■"(/t,/0=l-/^'(A,e) (2.6) 


As the process is iterated, these probabilities must stay in the range from 0 to 1. The 
formulas given above are the basic expression of the relaxation technique. This concept 
can be expanded and modified for other applications. This technique has the following 
characteristics: 

1. High compatibility pixels tend to reinforce each other and the low compatibility 
pixels tend to discourage each other. 

2. The degree of reinforcing oi discouraging done by a neighbor is proportional to its 
own probability of assignment of that class. 

3. The probabilities p{k,s) remain in the range from o to 1 and sum to 1 over all of 
class s. 

B. APPLICATIONS OF RELAXATION METHOD 
I. Line enliancenien t 

One of the applications in which relaxation techniques have been proven useful 
is the detection of long smooth curves in an image. Initially, line, detection operators 
are applied to the image, and their outputs are used to determine an initial probability 
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of the point lying on a curve with a given orientation. Tliese probabilities are then iter¬ 
atively reinforced: the probability of a point lying on a curve is reinforced by the prob¬ 
abilities of the other points lying on curves that smoothly continue it. After iterations 
of the procedure, points that lie on smooth curves tend to have high probabilities of 
being curve points, while the other points do not. 

An application of the relaxation technique to the detection of major edges in 
an image will be presented in the next chapter in detail. Initially, a gradient edge opera¬ 
tor is applied to the image. This provides information about edge strength and orien¬ 
tation at each point. These probabilities arc then iteratively reinforced. 

The general idea of the edge reinforcement process is as follows. 

a. Magnitude and direction 

First, the magnitude and direction of the image gradient are computed for 
each pixel k. The magnitude of pixel /c, divided by the maximum of the magnitudes over 
the entire image, defines the initial probability p%k,s) of a pixel k being an edge. 

b. Reinforcement of edges 

The reinforcement process defines a new edge probability of a pixel k in 
terms of the old probabilities at the pixel k and its ncignbor pixels. The computation of 
the new edge probability can be broken into steps as follows. 

(]} Interactions between edge and non-edge. The interaction between 
edge and edge, edge and non-edge, non-edge and edge, non-edge and non-edge depend 
on the edge direction at the two positions, on the orientations of the line joining the 
points, and on the distance between them. 

(2) Computation of the new probability. To compute the new edge prob¬ 
ability obtained from the iterations, first, compute weighted sums of the edge and non- 
edge probability increments, q'{k,e) and q'{k,n). These sums .arc then normalized to lie 
in the rangc(-l,+ 1) according to (2.4). They are used to update the edge and non-edge 
probabilities, while the updated values arc normalized so that they can sum to 1. 

The real issue in this thesis is to search for good line detection tech¬ 
niques applicable to lolargrams,and not to edge detection. Therefore, we are going to 
concentrate on the reinforcement of the relaxation method as applied to lofargratn. 

2. Labeling 

There is another application using the relaxation technique, which is called "la¬ 
beling." The objective is to assign labels to objects in an image. [Ref. 10] 
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a. Components of labeling 

In a finite set of relations between objects, the objects usually correspond 
to entities to be labeled. The objects are often geometrically or topologically related to 
each other. An input scene is thus a relational structure of all objects. In the simplest 
case, each object is to be assigned with a single label. Labels may be weighted with 
"probabilities" indicating the "probability of an object having that label. "Constraints 
determine what labels may be assigned to an object. A basic labeling problem is then: 
Given a finite input scene (relational structure of objects), a set of labels, and a set of 
constraints, find a consistent labeling. 

There are several types of labeling using relaxation. 

b. Discrete labeling 

It is a parallel iterative algorithm adjusting all object labels. All possible 
labels to each object arc assigned in accordance with constraints. Iterations are per¬ 
formed until a globally consistent labeling is found. In parallel, from each object's label 
set all the labels that are inconsistent with the current labels of the rest of the relational 
stiuctuie are climinatecl. 

c. Linear labeling 

The labeling process starts with an initial assignment of weights to all ob¬ 
jects. Weight are reminiscent of probabilities, reflecting the "probability that a label is 
correct." ilere p refers to probabilily-like weight rather than to the value of a probability 
density function. Let a relational structure with n objects be given by <?*, /c = 1 ,..m, each 
with nt discrete labels and let p^{?.) denote the weight, or the "probability" that the 
label is correct for the object a* and restiictcd as follows 

0 <:/>,(;.) ^1 ( 2 . 10 ) 

DrtW-i Pii) 

.1 

The linear labeling operator is based on the compatibilities of the labels, 
which serve as the constraints, A compatibility p^, looks like a conditional probability. 

( 2 . 12 ) 

A 
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The pyf^k I >T) may be interpreted as the conditional probability that object o* has label 
X given that another object a, has label X'. The operator iteratively adjusts label weights 
in accordance with other weights and the compatibilities. A new weight p^iX.) is com¬ 
puted from old weights and compatibilities as follows. 

/ 

The c„ are coefficient such that 




d. Nonlinear labeling 

Ii the compatibilities are allowed to take on both positive and negative 
values, we can express strong incompatibility and strong compatibility. Denote the 
compatibility of the event "label X on aC with event "label X' on a,” by >!')• 
two events occur together ol'ten, r*, should be positive. If they occur together rarely, 
should be negative. If they are independent, r*, should be 0. ’I'he compatibilities are based 
on correlations 



(2.15) 

a'(A') = /;(A-,A-)-(p(A0)' 

(2.16) 


(2.17) 


The r^| is used to obtain the positive or negative change in weight. 

I A- 

The q\{X) is the increment of the weight. The probability change is as follows 
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(2.19) 


Basically, labeling is diUbrcnt with the detection problem. Labeling is a 
process to assign the correct label to the correct position while detection is a process to 
take the desired portion of a scene having a mixed feature. Therefore, labeling is a dif¬ 
ferent technique for line detection. 

The general concept of the relaxation technique and its applications were 
discussed in this chapter. In the next chapter, line and edge detection algorithms using 
relaxation technique will be discussed in details. 
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III. RELAXATION ALGORITHMS APPLIED TO LINE DETECTION 


There are many algorithms already developed for line detection, in this thesis, the 
objective ij to select and test some of them for the lofargram processing. 

A. EDGE ENHANCEMENT(RXEG) 

This algorithm is derived from the original idea developed by Schachter, et al to re¬ 
inforce the contin.mus edges detected by a gradient operation. RXEG is an acronymn 
for "Relaxation method for EdGe dctcction".Edges reinforce other edges that interact 
with nearby non-edge points in specified ways. The gradient of the edges are also iter¬ 
atively adjusted. (Ref. 11) 

1. Initial edge probability 

A digital gradient operation is applied to the given image f. if we denote the .t 
and y components of the gradient by A/and A/, then the magnitude and direction of 
the gradient .ire given by 

(3.1) 


0 = tair'( 



(3.2) 


'1 he tmg of die gradient indicates the strength of the edge. 1 he angle 0 of the gradient 
is perpendicular to the edge direction. Edge direction is obtained by adding 90' to the 0 
angle of the gradient. 

We dclined the "probability" of an edge at a given point k by 


p{k,(!) 


mag{k) 


(3.3) 


where the max is taken over the entire image. The probability of a nonedge at a pi.xel k 
is defined as p{k,n) — \—p[k,e). 

2. Pixel and neighbor interaclion 

The centered pixel k could have four kinds of interactions with its neighbor pixel 
/ such as edge/edge, edge/non-edge, ncn-cdge/cdgc and non-edge/non-edge. 


13 





0. Edgejcdge inteiaction 

k is the centered pixel position, and / is one of the neighbor pixel position. 
Let a be the edge direction at k, P the edge direction at /, y be the direction of the line 
joining A to /, D the chessboard distance from k to /, i.e, max( 1 k^-l, |, | l<y~ly 1). Then 
the edge/edge reinforcement process between the points k and / has strength given by 

c{k,e\l,e) = cos(a-y) cos(/?-y)/2” (3.4) 

To sec the significance of this definition, a few simple examples are consid¬ 
ered. In these examples, the arrows indicate the direction along the edge, with the dark 
side of the edge on the left. These examples showed that parallel and perpendicular 
edges have no effect on one another. 


Case 

a 

P 

y 

cos(a - i) cos(p - Y )/ 2 ‘^ 

n 

90 

90 

0 

0 


90 

270 

0 

0 

A 

A 

90 

90 

90 

1/2^ 

A 

90 

270 

90 

-1/ 2^ 


90 

0 

0 

0 

■■ 

90 

180 

0 

0 


Table 1. Examples of edge / edge interaction [Ref. 8] 


b. Edgejnon-edge interaction 

Besides the edge/edge interaction which occurs between the edge probabili¬ 
ties p{k,e) and j){l,e) there are also interactions involving the non-edge probabilities 
p{l,n). The edge probability at pixel k is weakened by the non-edge probability at pixel 
/ to the degree c{k,e\l,n) defined by 

c{k,e-,l,n) = min[0, - cos(2a-2y)/2^] (3.5) 
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In the examples, non-edge neighboring points collincar with the center edge points 
weaken them, whereas non-edge points alongside edge points have no eflect on them. 


Case 

a 

y 

- cos(2a “ 2 y) 12^ 

c(k, e; 1, n) 


90 

0 

H2^ 

0 


270 

0 

112^ 

0 

• 

90 

270 

-\I2^ 

-1/2*^ 

A 

90 

90 

-H2^ 

-112^ 


Tabic 2. Examples of edge / non-edge interaction (Ref. 8] 


c. Non-edge I edge interaction 

'1 he non-edge probability at pixel k is alTectcd by the edge probability at the 
neighbor pixel / to the degree c[k,n\i,e) defined by 


c{k,n\i,e) - 


(1- cos(2/?-2r)) 


(3.6) 


Case 

a 

P 

y 

(1 -cos(2a - 2y)) /2^ 

■A 

0 

90 

0 

1/2° 

A 

0 

90 

90 

0 


Tabic 3. Examples of non-edge / edge interaction (Ref. 8] 


I'rom the examples, neighbor edge points alongside the center point strengthen them, 
While edge points collincar with non-edge points have no eflect on them. 
d. Non-edge}non-edge interaction 

Ihe non-edge probabilities at centered pixel k and the neighbor pixel I re¬ 
inforce each other to the degree c{k,n',l,n) defined by 

c{k,ir,l,n) = — (3.7) 

2 
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Non-edge probabilities arc reinforced by the otlier nearby non-edge probabilities. 
c(/-,«;/,«) is directionally independent. 

3. Combined reinforcement proce.ss 

For each pixel k, the net elfcct of its neighboring pi,xcl on its edge probability 
p{k,e) and non-edge probability pik,n) = is computed as follows 


q'{k,e) = 

/ 


+y ^C2p'il,it)c[k,c;l,it) 


(3.8) 


= Y^Qp%<!)c{k,ivJ,e) 

I 


I 


where C,, C\, C], Q are constants whose sum is equal one. The standaid values of 
C, = 0.86(5 , Cj = O.I24, 0, = 0.005, and C« = 0.005 was used in the paper by Bruce 
J.Schaciiter, et al [Ref. 11]. The results of the iteration process are somewhat sensitive 
to the choice of the Cs. For example, if C, is too large, the edge will thicken and will 
be extended into non-edge points; while if C 4 is too large, gaps will appear at W'eak spots 
in the edges and at sharp angles. 


q'(k,e)= -r - 

^ I q{k,c) I -t-1 q\k,n) I 

- 

I q\Ke) I -1-1 q\k,n) \ 
p'{k,e) =p\k,e)l\-\-q'{k,e)'] 
p'{k,n)== p\k,ii)iU-q\k,n)] 


(3.10) 

(3.11) 

(3.12) 

(3.13) 
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(3.14) 




p'{k,e)+p’(k,it) 


This process is then iterated by using p'^'{k,e) iti place of p'{k,e), and by using 
in place of p'{k,n). Wc also compute the estimated edge direction A',(A) and 
at each point 

A'#)={iy(/.c)cos(0(/<)) 


+^p\l,e)c\k,c\l,e) cos(0(/)) 


/ 


A'^,(/c) = ]Vp\l,e) %m[0{k)) 


(3.15) 


+^//(/.c)c(/v,<?;/,r) sin(0(0) (3.16) 

I 

r'(A) = tair'(A',(.v^')/A>j-)) (3.17) 

For large values of the constant W, is close to 0; while for small values, it is strongly 
inllucnced by the neighboring O's. 

n. LINE AND CURVE ENIIANCEMENTlRXLN) 

In section A, only two categories arc associated with the image pixel, (i.e, edge and 
no.i-cdge). In this section, more categories arc considered in the RXLN relaxation pro¬ 
cedures. RXLN is an acronymn for "Relaxation method for LiNc detection." 

This algorithm was devlopcd by Zuckcr, Ilununcl and Roscnfeld (Ref. I2J. The re- 
la.xation process i> applied to the detection of .smooth lines and curves in noisy images. 
Nine class labels associated with each image point will be considered. Light classes in¬ 
dicate lines at various orientations and one indicates the no-Iine case. In the relaxation 
process, interaction takes place between the probabilities at neighboring points. This 
permits line segment with compatible orientations to strengthen one another. Similarily, 
no-line class is reinforced by neighboring no-line class and weakened by the neighboring 
oriented line classes. 
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1. Initial piobabilities 

At each point, probabilities are assigned to the nine classes which are lines in 
eight possible directions and a no-linc class. The initial probability for each class is ob¬ 
tained by evaluatng a nonlinear line detector as shown in Figure 6 at every picture point 
along the eight orientations. Eight sets of pixels representing eight orientations are se¬ 
lected and shown in Figure 6. Each set of pixels arc rearranged as the template shown 
in Figure 6 (b). When the condition in Figure 6 do not hold, the response is zero. When 
they do hold, the response is calculated by 

R=={D+E+H)-~{A+D+G+C+F+r) (3.18) 

At each pixel, the detector's response are computed for every orientation, 'i hen only one 
orientation which has the maximum detector's response is selected as the orientation of 


the pixel. [Ref. 13] 






.Temalala 


Conditions 


A 

B C 


A 

< B > C 


D 

E F 


D 

< E > F 


G 

H 1 

(a) 

G 

< H > 1 


B B 

B 

B 


H 

H 

E E 

E 

HE H 

E B 

HE E 

E 

H H 

H 



B B 

B 

(1) (2) 

(3) 

(4) 

(b) 

(5) 

(6) (7) 

(8) 


Figure 6. Nonlinear line detector, (a) Nonlinear line dctecioi for the vertical ori¬ 
entation. When the conditions do not hold, the response is zero; when 
they do hold, the response is calculated by (3.18). (b) Eight orientations 
of the detector. 
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To describe this process, let d,(k) denote the output of the nonlinear detector for 
orientations s at pixel k . Let p{k,s) denote the probability that a pixel k belongs to class 
s. Then the initial probabilities can be obtained by scaling the detector's responses at 
each position by the maximum over tlic entire picture. The probability that pixel k has 
class s is expressed as follows 


p{k,s) 


nvxMk)) 

j=i _ 

max((max)(/j(A)) 

S=\ It 


8 

1=1 


(3.19) 


max(rf,(^)) means the one selected maximum value over all orientations. max(<y,(/0) 

j-i * 

means the maximum value of oriented s over all image pixels, and the probability that 
pixel k has no-line is defined by 


8 


p{k,9)= ]-^_j{k,s) 


1=1 


(3.20) 


2. Conipalibility coefficients 

To update the probabilities, the compatibility relation between a pixel and its 
neighbors must be specified.. The compatibilities depend only on the relative orientations 
of the neighboring point. These compatibilities can be specified in the following way. If 
two neighboring line segments arc oriented in the same direction or close to the same 
direction, they add support to one another. If, on the other hand, two neighboring seg¬ 
ments arc oriented perpendicularly to one another, they will reduce support to each 
other. All other pairs of line segments are distributed between these two extremes. In this 
algorithm, there arc three types of compatibility coefficients. 
a. Compatibility coefficients between lines 

Figure 7 shows the compatibility coefficients between the lines in five ge¬ 
ometrical relationships. As standard values, from left to right: 1.0, 0.5, 0.05, -0.15, -0.25 
are used. When straight line enhancement is desired, 1.0, 0.0, -0.1, -0.17, -0.22 ate used. 
[Ref. 13] 
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b. Compatibility coefficients between fine and no-lines 

As its standard value, -0.1 is used. Decreasing the value to -0.5 increases the 
elTect to fill the gap between line segments. [Ref. 13] The standard value is 0.25. 

I / / 

'■0 .5 ,05 -.15 -.25 

Figure 7. Compatibility weights between line labels 

f. Compatibility coefficients between no-lines 

Increasing the value up to 0.5 improve the noise efiect but, may cause era¬ 
sure of line labels. [Ref 12] 

For each pixel /<, its neighbor pixel has 9x9 cases of interaction with the 
center pixel. If neighborhood size is 5.x5. ’Ihcrc are 5x5x9,x9 cases of compatibility coef¬ 
ficients, 

c(/£,l;l,l) . . c(A-,l;25,l) 

* • • • 

• • • • 

c(/<,9;l,9) . . c(*,9;25,9) 

3. Updating probability 

The updating proccs.s can be expressed in terms of compatibility functions. Let 
q'{k,s) represents the increment applied to p'{k,s), which is the probability that pixel k 
belongs to class s. The probabilities of each pixel in each class are computed as follows 

= Yj ( 3 - 21 ) 

i t 



where s repr^'sents nine classes of reference pixel and t represents nine classes of neigh¬ 
boring pixels and 









(3.22) 


r+V/ s P%s)li-\q\k,s)'] 

^p\k,s)[l+q\k,s)^ 

S 

C. LINE AND CURVE ENHANCEMENT(RXLA1) 

This algorithm is delivered from the original idea developed by S.Peleg and 
A.Rosenfeld (Ref. 14). RXLAl is an acronynin for "Relaxation method for Line And 
curve." Basic concept using in this algorithm is the same as that used in the RXLN(sce 
section B) except that the method for compatibility cociricient is different. 

1. Compatibility coefficient 

One possible interpretation of the compatibilics is in terms of statistical corre¬ 
lation. Couclation has properties as follows 

a. If class s and class t arc compatible for pixel k and pixel /, c{k,s\l,t) > 0 

b. If class s and class t arc incompatible for pixel k and pixel /, c{k,s\l,t) < 0 

c. If neither class is constrained by the other, c{k.s',t,t) = 0 

d. The magnitude of c{k,s\l,i) repicsents the strength of the compatibility. 

Estimates of the correlation cocHicicnts derived from analyzing the initial 
probability are 


c{k,s\l,t) 


a(s)aii) 


(3.23) 


/'(/{,s) is the initial estimate of the probability of pixel k with class s, p{k,s) is the average 
of p(/<,s) for all pixels k, a(.s) is the standard deviation of p{k,s). If the neighborhood size 
is 5 X 5, 5 X 5 X 9 X 9 cases exist for each pixel as below 

c(/<,I;l,l) . . c(/v,l;25,l) 

• » • • 

• • • • 

c(^,9;l,9) . . c(*,9;25,9) 

Then, for every case and for all pixel /< we have the summation as given by 
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c(A,l;I,l) . . c(/<,l;25,l) 

• • • • 


c(A„l;l,l) . . c(/<„l,25,l) 

• • * * 

• • • • 

C(A,9;1,9) . . c(/<,9;25,9) 


• • • « 

c(A„9;l,9) . . c(A„9;25,9) 


c(A-2,1;1,1) • . c(A2,1;25,1) 


+ 


c(/<2.9;1,9) . . c(/<2.9;25,9) 


+ • • • 


Thcrerorc 




k _ _ 


a(5)a(f) 


(3.24) 


In tliis algorithm, there is the cllect of dominance, that is, wlien one class 
dominates the picture, its correlation coellicicnts with all other classes arc high. Thus in 
the relaxation updating, the dominant class gets most of the support, and after a few 
iterations almost all the pixels have this class. '1 he clfcct of dominance among classes 
can be alleviated by weighting the c(A,.v;/,/) by the probabilities that the corresponding 
classes do not occur. This greatly reduces the values of the coellicicnts involving domi¬ 
nant classes, but will only slightly reduce the values of the coellicicnts involving rare 
classes. '1 he compatibility coellicicnts used in this algoritm arc 

c(A,s;/,0 = C1 -p(/<,5)][1-/T(/,/)] X c{k,l,s,t) (3.25) 


Two line and edge detection algorithms using relaxation technique were dis¬ 
cussed in this chapter. These algorithms will be applied to the lofargram for line de¬ 
tection. Experimental results of these algorithms will be discussed in Chapter IV. 
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IV. EXPERIMENTAL RESULTS AND A MODIFIED RELAXATION 

APPROACH 


In the experiments, four difFerent relaxation algorithms are applied to lofargrams 
having difFerent S/N. Both real and artificially created lofargrams were used as test im¬ 
ages. As discussed previously, two algorithms RXLN and RXLAl arc line detection al¬ 
gorithms while the third algorithm RXI2G is an edge detection algorithm. 1 he purpose 
of the experiment is to test these three algorithms to sec how they work for lofargrams. 
■fhe secondary purpose of these experiments is to obtain thin lines for loFargram inter¬ 
pretation. Because the results of these algorithms were not satisfactory', improvement 
oFthc most suitable algorithms arc suggested. 

The results of the line detection algotithms RXLN and RXLAl were not encour¬ 
aging when the S,'N is low. 'Fhe reason will be explained later in this chapter. Fhe results 
improve as the S/N increases. One effect is that the noise was removed very slowly in 
iterations of reinforcement. On the other hand, the edge detection algorithm, RXHG 
showed better results than the above two algorithms in the presence of low S/N. But the 
results were still not ideal. A double line was shown when the S/N was low’ because 
RXLG was the edge detection algorithm. 'Fheieforc, some modification is needed for line 
detection in tlie RXIiG. This will be di.scusscd in details later. 

Two parts in the RXIIG algorithm were modified. The first one is the initial magni¬ 
tude of the pixel, since the outputs arc greatly influenced by the initial magnitude of the 
input image. I hc other is the initial direction of the pixel. Because RXIIG uses a difl'cr- 
ential type of edge detection operator, the initial direction in RXEG is not appropriate 
for the lofargrarn line detection. In RXIIG the edge direction has a dark side which in¬ 
dicates the outside of the object and a blight side which indicates the inside of the object 
body. But the direction of a line simply indicates the potential direction of a track in a 
lofargrarn. Detail will be discussed in section C. 

In a search for a better initial magnitude and direction, many other algorithms w'ere 
tested. It assures that the magnitude and diicction output from EGPR which is another 
line detection can be useful. The RXEG was modified and tested accordingly. 

A. LINE DETECTION ALGORITHMS (RXLN, RXLAl) 

Both RXLN and RXLAl are line detection algorithms which use similar methods 
for initial probability and compatibility coefficient calculations as explained in Chapter 
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HI. But the applications of these are very discouraging. Figure 8 shows a real lofargram 
as the input image and its initial magnitude calculated in RXLN. Figure 9 shows the 
results of two iterations of reinforcement. And Figure 10 shows an artificial lofargram 
with S/N of 6 db and its initial magnitude image. Figure 11 show the results of the four 
iterations of reinforcement. And Figure 12 show the results for 3 db and 6 db of S/N. 

The initial magnitude of input image shows that the result is vciy noisy. In these 
algoritlims, there is the eficct of dominance. When one class dominates the picture, high 
value is assigned to the dominating class at each pixel and low values are assigned to the 
other classes. Thus, the dominant class gets the most support in the relaxation updating 
process. And after a few iterations almost all the pixels have the same class. Figure 8 
shows tliat the noise dominates the picture when the lofargram has low S/N. Figure 10 
(b) still shows that the noise dominates the picture. When the S/N is higher, the signal 
showed up weakly as shown in Figure 11. But noise still dominates the picture. Gener¬ 
ally noise dominates the lofargram. Noise sttongly dominates the lofargram when the 
S/N is low. Hence the RXLN and RXLAl algorithms are not suitable for lofargram 
processing. 

B. THE RESULTS OF THE EDGE DETECTION ALGORITHM RXEG 

Figure 13 shows the initial magnitude of the RXBG of the input image. Figure 14 
shows the results of the reinforcement applied to Figure 13. In this implementation, 
neighborhood size is 5x5, then the 24 neighboring pixels inlluenced the centered pixel. 
Figure 15 shows the artificial lofargram input image and its initial magnitude with 6db 
of S/N. The iteration scheme has considerable noise removal power as shown in the 
example of Figure 14. Figure 16 show's that the results of the artificial lofargram with 
higher S/N can yield a stronger response. A double line in Figure 14 (b) and Figure 16 
indicate that the edge w'as detected from the input image. Because RXFG uses the edge 
detection operator, double Hires are shown at lower iterations, The weaker line of the 
double lines disappeared gradually in the lelaxation process w’ith 4 successive iterations 
as shown in Figure 14. Figure 17 shows the results for 3db and 6db of S/N. A desirable 
example of edge detection is shown in Figure 18. 

C. MODIFIED EDGE DETECTION ALGORITHM RXEG 

The results of RXEG shows that tiiis algorithm works as an edge detection algo¬ 
rithm. It can detect tlie edges in the lofargram, even though the S/N is very low. Double 
lines in Figure 15 (b) shows an example of the edge detection. 
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Figure 8. 



(*) 



(b) 

A real lofargram input image and its initial magnitude from RXLN 
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Artificial lofargram of 6dh and its initial magnitude from RXLN 









































Figure 13. Initial magnitude of Fig 8 (a) in RXEG. The Prewitt operator is used 
for initial probability. 
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(b) 

Figure 14. Four iterations (1,2,3,4) of RXEG 
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(a) (b) 


Figure 15. Artificial lofargram Mith 6db S/N and its initial magnitude 



(c) (d) 

Figure 16. Results of four iterations(l.2,3,4) of reinforcement from Fig 15 
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Figure 17. Results of tuo tofargrams, (a) 3dl) (b) 6db. 


Results of the reinforceiiient with 5 iterations. However, the purpose of lofargram 
processing is to detect the line of the signal for the interpretation of the lofargram. The 
desirable ideal example of the line detection is ^hown in Figure 19. To accomplish this, 
some changes in RXHG are required. Many algorithms were tested to generate the edge 
magnitude and line direction for the initial probability. The new magnitude and direction 
obtained from EGPR were adopted in the modified algorithm. 



(a) (b) 

Figure 18. Example of edge detection process (a) input (b) output 


1. Edge magnitude output 

The gray level of the edge enhanced magnitude of each pixel are processed by 
masks. Nine masks are formed around a pixel in the shape of pentagon, hexagon and 
rectangular as shown in Figure 20. The average values and the variances within the 
masks around a pixel are calculated. 'I hen. the average value of the mask with minimum 
variance is assigned as the gray level of the center pixel. 



















(a) (b) 

Figure 19. Example of line detection process (a) input (b) output 
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Figure 20. Nine masks from EG PR 



2. Line direction output 

The initial line direction of each pixel is determined from the direction of the 
mask with the minimum variance. Eight line directions can be associ itcd with each im¬ 
age point. The directions are symmetric to the origin. Therefore, only four of them are 
considered as 0’ = 180’, 45’ = 225’, 90’ = 270’ and 135’ = 315’. The direction for an edge 
is the direction along the edge with the dark side of the edge on the left and the bright 
side on the right side. This means that the edge direction 0’ and 180’ arc different di- 
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rections. However, the directions for line detection from masks in EGPR do not have 
dark side and bright side. The direction in EGI’R simply implies the direction of a line 
without dark side and bright side. This means that the line direction 0’ and 180" are the 
same direction as shown in Figure 21. The edge direction with bright on the right side 
arc shown in Figure 22. 



Figure 21. An example of line direction 

A 

dark / bright 

/ 

bright f dark 

Y 

Figure 22. An example of edge direction 

I'he four directions generated from the EGI’R are shown in Figure 13. But, start with 
the second iteration after the first rela.xaiion process, infinite line directions are consid¬ 
ered at each pi.xel. 

3. Pixel and neiglibor interaction 

Edge detection in RXEG detects double edges of the signal line in the lofargram. 
It can detect the edges reliably when the S;N is high. But, the purpose here is to detect 
thin line as shown in Figure 19. For the thin line detection, it is necessary to construct 
a compatibility function with proper compatibility values for the interactions between 
different classes. Basically, the same compatibility function as that of the RXEG can be 
used for the interactions betw'een lines. The only difierence between the line direction 
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90 



0 

Figure 23. Four possible line direclious 

and the edge direction is that the directionality of a line greater than 180’ is considered 
as the same as its symmetric value througlt tite origin. In other words, the compatibilities 
of lines between tiie center pixel and the ncigltboring pixels can be handled in the same 
way as in RXIiG. 

'I he center pixel k could have lour kinds of interactions with its neighbor pixel 
/, i.c, linc/linc, linc/no-linc, no-linc/linc and no-linc/no-linc. 
rt. Linclliite interaction 

Let k be the centered pixel position, and / the neighbor position. Let a be 
the line direction at k, (i the line direction at the neighbor pixel /, y be the direction of 
line joining k io I, D the chessboard distance from k to /, i.c, ntax( | k-l, |, 1 k,-ly |) . 
Then the li le/line reinforcement process between k and / has strength given by 

c(/<,/;/,0 = cos(a-y) cos(/?-y)/2^ (4.1) 

which is similar to eq (3.4). Here / represents the line class and n represents non-line 
class, 1 0 sec the significance of this definition, a few simple examples are considered in 
Table 4. Results changed with did'erent compatibility value of parallel lines. In the ex¬ 
periment, when a higher value of compatibility is chosen for parallel lines, background 
noise increases, while a lower value is chosen when the opposite situation takes place. 
b. Linclno-line interaction 

'I he line probability at pixel k is weakened by the no-linc probability at pixel 
/ to the degree c{k,l’,l,n) defined by 

c{k,l]l,n) = minCO, - cos(2a-2y)/2^] (4.2) 

which is similar to eq (3.5). A few' examples are shown in Table 5. 
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c. No-UnelHm interaction 

I'he no-iinc probability at pixel k is alTectcd by the line probability at the 
neighbor pixel / to the degree c{k,n\l,t) denned by 


c{k,iv,l,l) = 


(l--cos(2/?-2v)) 


>n+i 


(4.3) 


which is similar to eq (3.6). A few examples are shown in Table 6. 


case 

a 

P 

7 

cos{a - y) cos(P -y) 

1 

1 

90 

90 

90 

1 / 2 '^ 

■Hi 

90 

0 

0 

0 

n 

90 

90 

0 

— 


Tabic 4. Examples of line / line interaction 


case 

a 

Y 

- cos(2a - 2y) / 2 *^ 

c{k, 1; 1, n) 

n 

90 

0 

\I2^ 

0 

n 

90 

90 

- 1 / 2 ’^ 

- 1 / 2 '^ 

n 

90 

270 

- 1 / 2 ” , 

-1 /2’^ 


Table 5. Examples of line / no-line interaction 


case 

a 

P 

Y 

(1 - COS(2P - 2y)) / 2 '^ 

•1 

0 

90 

0 

1 / 2 '^ 

t 

1 

0 

90 

90 

0 


Table 6. Examples of no-Iine / line interaction 
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d. No-lmcino-linc interaction 

The no-Hne probabilities at the centered pixel k and that of the neighbor 
pixel / reinforce each other to the degree c{k,n\l,n) defined by 

c{k,n\l,n) = -\ (4.4) 


which is similar to eq (3.7) 

4. Combined leinfui cement process 

For each pixel k, the net effect of its neighboring pixei.. an its probability p{k,[) 
and no-linc probability p(/<,«) = \-p{k,l) is computed as follows 


I 

/ 


(4.5) 


qih.n) = 

/ 


+^Qp'{l,n)c{k,n;l,n) (4.6) 

I 


where C,, C„ Q, Ci are constants who.sc sum is taken to be one. The standard values 
needed here are C, = 0.866, Cj = 0.124, C, = 0.005 , and Q = 0.005. 'fhe results of the 
iteration process are somewhat sensitive to the choice of the C's. For examjile, if C, is 
too large, the line will thicken and will be extended into no-linc points; while if Q is too 
large, gaps will appear at weak spots in the lines and at sharp angles., 


qW) = 


I I +1 I 


(4.7) 


q'{k,n) 


I 7V<.0 I + I (i\k,n) I 


(4.8) 
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p\k,[)=p\k,l)U+q'{k,l)-] 

(4.9) 

p'{k,n)=p\k,n)l\+q'{k,n)'] 

(4.10) 

P p\kMp\k,n) 

(4.11) 


This process is then iterated with replacing and l-p'*'{l<,f) replacing 

p'{li,n). We also compute the estimated line direction at each point 

cos(0(/()) 

+Y^p\lMkAl,f)cosm) (4.12) 

i 

A'y{k) - U'p\k,l) m\{0{k)) 

+^/>V/)c(A,/;/,/) sin(0(0) (4.13) 

I 

0^+'(/<) = tan-'(A',(/<)/A',(/<)) (4.14) 

where IV is a constant. For large values of the constant W, 6'*' is close to 0; while for 
small values, it is strongly influenced by the neighboring 0's. 

5. Experimental results 

Figure 24 shows the initial magnitude of the input image of Figure 8 (a). Figure 
25 shows the results for 3db and 6db of S/N. I'igure 26 shows the results of four iter¬ 
ations of the reinforcement applied to Figure 24 using the modilied RXEG. Figure 27 
shows the artificial lofargram with 6db of S/N and its initial magnitude. Figure 28 shows 
the results of four iterations of reinforcement with 6db S/N artificial lofargram. Com¬ 
paring Figure 13 and Figure 24 reveals that the initial magnitudes of the modified RXliG 
are improved by using the EGI’R line detection. The initial magnitude of modified 
RXFG shows stronger signal line in Figure 24 when compared to Figure 13. However, 
the end results of reinforcement do not improve significantly. From Figure 28, it is ob¬ 
vious that a double line still exists and that the lines are thicker than those of the RXEG 







algorithm. In this experiment, the purpose is to obtain single thin line from the 
lofargram. But the result is not fully satisfactory. A few things could be dis issed with 
respect to the suspected reasons of the double lines and of the thicker line in the result. 
The initial magnitude as shown in Figure 27 (b) was caused by the Prewitt operator of 
the RXEG. In Figure 1.5 the result of RXEG thows strongly the splitting line. Therefore, 
the reinforcement process couldn't be the problem of causing double lines. Figure 29 
shows the diagram of the modified RXEG. As for the possible cause of the thicker line, 
the modified initial magnitude may be the reasot.. Comparison of Fig 15 (b) and Fig 27 
(b) reveals the improvement of the initial magnitude. The initial magnitude of the mod¬ 
ified RXEG shows stronger and thicker lines. Consequently, the results of the re¬ 
inforcement at a later stage applied to this initial magnitude shov;*' thicker lines. 
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Figure 24. Inilial magnitude of Fig 15 (a) In modified RXEG 


To solve the double line problem, it is necessary to find a method to yield the 
initial magnitude in single lines. Two dilfeicnt method were tested for this purpose. In 
the first method each pixel was divided by the maximum value over the entire picture 
as follow. 
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(4.15) 


rikJn) 

max p{k,ln) 

k 

This method removed noise very quickly. The initial magnitudes at each pixel are cal¬ 
culated by dividing each pixel by the maximum values over the entire picture. Therefore, 
very small values are assigned to noise and the high values are assigned to the signal line 
when the S/N was high. On the other hand, when S/N is low, the high values are as¬ 
signed to both the noise and the signal because the difference of gray level between the 
noise and the signal is small. As a result, this method is still not suitable to lofargram 
processing. In the second method, each pixel was divided by the maximum value of the 
local neighborhood as follow 


pjkjn) 

maxp{ljn) 


(4.16) 


Initial magnitude in this method depends on the local maximum value. In the exper¬ 
iments it is realized the initial magnitude could be white or dark entirely depending on 
the characteristic of the noise of the lofargram. This second method was not appropriate 
for the calculation of the initial magnitude before the relaxation process. 

In this chapter various results of the algorithm studied were discussed. Thick 
line problem were discovered. Several attempts has been made to solve this problem. 
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(d) 

Fioure 26. Four iterations (1,2,3,4) of modified RXEG (continue) 
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(a) (b) 

Figure 27. Arlificial lofargram >^ith 6tib S/N and its initial magnitude 



(c) (d) 


Figure 28. Results of four iterations (1,2,3,4) from Fig 27 in modified RXEG 




























Figure 29. The procedure of niudified RXEG 
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V. CONCLUSIONS AND RECOMMENDATION 


The fact that each algorithm has its own peculiar characteristics in the lofargrani 
processing was found in this study. 1 hese characteristics arc very much dependent on the 
S/N of the input lofargrani. Algorithms RXLN and RXLAl using similar ideas can de¬ 
tect lines in the loHirgram when the S/N is relatively high. But, the result were very dis¬ 
couraging when the S/N was low and noise smoothing was done slowly in the iterations 
of the reinforcement. Our main concern of the lofargrani processing is to detect a thin 
line when S/N is low. ‘Ihcrcforc, thc.se algorithms arc not suitable for lofargrani process. 
On the other hand, the RXilG algorithm results in good lofargram processing. The only 
shott coming is that it is an edge detection algorithm. It smoothes the noise and en¬ 
hances the line components in the lofaigram even at low S/N. The last algorithm 
studied in this thesis is the inodincd RXIid algorithm. Several method were tested to 
solve the double line problem. But, double line problem still e,\ists in the modified algo¬ 
rithm. ‘1 0 solve this problem, it is necessary to search for an algorithm that can yield thin 
single line in the image of the initial magnitude. 
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APPENDIX. THE PROGRAM OF MODIFIED RXEG 

INTEX3ERM ISX/256/, ISY/256/, ITER/4/, ISNX/5/,ISNy/5/, lSW/1/ 
INTBGERM IP (256,256), IPRB( 256,256), JP (256,256) 

DIMENSION C(4),CTYPE{4,2),THET0(256,256> 

DIMENSION PRBEl(256,256),PRBE2(256,256),0(8) 

DIMENSION THETl(256,256),THET2(256,256) 

BYTE G(256) 

REAL LENGTH, H/3.0/,MAX,MIN 

OPEN(UNIT-1,NAME-'BEAMOO.BIN',TYPE-'OLD',ACCESS-'DIRECT', 

2 RECORDSIZE-64,MAXR£C-256) 

OPEN(UNIT-2,NAME-*TLOFAR20-4X5.DAT‘ ,TYPE-'NEW' ,ACCESS-'DIRECT’ , 
2 RECORDSIZE-64,MAXREC-256) 


DO 10 1-1,256 
READ(1'I)G 
DO 20 J-1,256 

IP(I,J)-G(J) 

1F(IP(I,J).LT.0) IP(I,J>-IP{l,J)+256 
20 CONTINUE 

10 CONTINUE 

CALL RXEG (IP, THETl, THET2, PRBEl, PRBE2,3SX, ISY, ISNX, ISNY, 1 TER, 
& ISM,TEMP,JP,THET0) 

MAX-PRBE2(1,1) 

MIN-PRBE2(1,1) 

DO 30 1-1,256 
DO 40 J-1,256 

IF(MAX.LT.PRBE2(I,J)) MAX-PRBE2(1,J> 

IF(MIN,GT.PRBE2(1,J)) MIN-PRBE2(1,J) 

40 CONTINUE 

30 CONTINUE 

LENGTH-(MAX-MIN) 

DO 50 1-1,256 
DO 60 J-1,256 

IPRB(I,J)-JNINT(((PRSF2(I, J ) -MIN)/LENGTH)* 24 3.) 

60 CONTINUE 

50 CONTINUE 

DO 70 1-1,256 
DO 80 J-1,256 

IF(IPRB<I,J),GT.127) THEN 
G(J)-IPRB(I,J)-25e 
ELSE 

G(J)-IPRB(I,J) 

ENDIF 

80 CONTINUE 

WRITE(2'I)G 
70 CONTINUE 

END 

SUBROUTINE RXEG(IP,THETl,THET2,PRBEl,PRBE2,ISX,ISY 
S ,ISNX,ISNY,ITER,ISW,TEMP,JP,THET0) 

C 

C 

CP Edge reinforcement by releucation method. 
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c 

CS CALL RXBG(IP,THBT1,THET2,PRBE1,PRBE2,ISX,ISY 

CS /1SNX,ISNY,ITER,ISW,TEMP) 

C 

CK RELAXATION, EDGE 


IP(ISX,ISY) ! 

Input image 

(I) 

THBT1(ISX,ISY)! 

Gradient direction after ITER-1 
iterations 

(W) 

THi:,T2(ISX,ISY)j 

Gradient direction after ITER iterations 

(0) 

PRBE1(ISX,ISY): 

Edge probability for each pixel after 
lTER-1 iterations 

(W) 

PRBE2(ISX,ISY)8 

Edge pr^>ability for each pixel after 
ITER-1 iterations 

(0) 

ISNX,ISNY ; 

Size of neighborhood being considered 
(odd number, for example: 5,5) 

(I) 

ITER : 

Number of iterations 

(I) 

ISW t 

Switch (1 - standard, 2 - noise reduction 
type 

(I) 

Reference 




(31! B. J. Schachter, A. Lev, S W. Zuc)cer, and A. Rosenfeld, "An 
a ppl ication of relaxation methods to edge reinforcement," 

IEEE Trans. Syst., Man, Cybern., vol. SMC-7, pp. 813-816, 

Nov. 1977. 

JULY 1979. PROGRAMMED BY K.SAKAUE 

(1) When computing the Initial values for edge probabifty and 
edge direction using an operator other than the Prewitt 
operator, rewrite lower-level subroutine RXEP. Any operators 
that ca obtain a differential (Dx,Dy) can be used. 

(2) The processing characteristic largely depends upon selection 
of coefficients Cl, through C4 (See the reference.) The two 
standard )clngs of values have already been set in array CTYPE 
us.'.ng data statements, and values can be changed with ISW. 
vr.' n using other values or changing values depending upon 
ti.t number of iteration tiroes, rewrite the program so a proper 
value is set to the argument C in lower-level subroutine RXEI. 
The following condition however, must be met. 

C(l)+C(2)+C(3)+C(4)-1 

(3) The coefficient W used when updating (THETl) is set in data 
statement to 3.0. The same value was also used in the ref feme 
This value may be changed as neccessary. 

DIMENSION IP(ISX,ISY),JP(ISX,ISY) 

DIMENSION THETl(ISX,ISY),THETO(ISX,ISY) 

DIMENSION THET2(ISX,ISY) 

DIMENSION PRBE1(ISX,ISY) 

DIMENSION PRBE2(ISX,ISY),TEMP(256,256) 

DIMENSION CTYPE(4,2) 

DIMENSION C(4) 

DATA CTYPE/.866,.124,.005,.005 

4 ,.706,.176,.059,.059/ 

DATA W/3.0/ 

CF INITIAL EDGE VALUES AND INITIAL PROBABILITIES 
C 
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CALL BGPR(IP,JP,1SX,ISY,THET1,THET0) 

CALL RXro(lP,THETl,PRBEl,lSX,lSY,TEMP) 

C 

CP RELAXATION UPDATING PROCESS 
C 

AMN-0. 

AMX-0. 

DO 9 1-1,4 

C(I)-CTYPE(I,ISH) 

9 CONTINUE 
DO 10 I-1,ITER 

IF(I.BQ.l)GO TO 11 
DO 12 IY-1,1SY 
DO 12 1X-1,1SX 

THETl<IX,lY)-THET2(IX,lY) 

PRBEl(IX,lY)-PRBE2(IX,lY) 

12 CONTINUE 
11 CONTINUE 

CALL RXEI (THETl, THET2, PRBEl, PRBE2, ISX, ISY, ISNX, ISNY,C, W) 
10 CONTINUE 
RETURN 
END 


C 

C 

CS 

C 

CP 

C 

CK 

C 

CA 

CA 


SUBROUTINE EGPR(IP,JP,ISX,ISY,THETl,THETO) 
Copyright (c) 1983 by AIST MITI(JAPAN) 

CALL EGPR(IP,JP,1SX,ISY) 

Edge preserving smoothing operation. 

EDGE PRESERVING SHOOTHING, SMOOTHING 

IP(ISX/ISY) : Input image array 

JP(ISX,ISY) : Output image array 


(I> 

( 0 ) 


M 

CM 

CM 

C 


C 


The program is considerably long, but the main part consists 
of about 120 lines in the first half and the remainder is used 
as routines for edge processing. 

DIMENSION IP(ISX,ISY),THETl(ISX,ISY),THETO(ISX, ISY) 

DIMENSION JP(I3X,ISY) 

DIMENSION A(9), V(9),0(9) 

DIMENSION K(5,5), L(5,5) 


EQUIVALENCE 

& 

& 

EQUIVALENCE 

& 

& 


(A(1),A1), (A(2),A2), (A(3),A3), 
(A(4),A4), (A(5),A5), (A(6),A6), 
(A(7),A7), (A(8),A8), (A(9),A9) 
(V(1),V1), (V(2),V2), (V(3),V3), 
(V(4),V4), (V(5),V5>, (V(6),V6), 
<V(7},V7), (V{8),V8), (V(9),V9) 


C 


EQUIVALENCE 

& {K11,K(1,1)),{K21,K(2,1)),(K31,K(3,1)),(K41,K(4,1)),(K51,K(5,1)), 





t (K14,K(1,4)> 
t (K15,K{1,5)) 
EQUIVALENCE 
& (L11,L(1,1)) 
(L12,L(1,2)) 
(L13,L(1,3)) 
(L14,L(1,4)) 
{L15,L(1,5)) 


/(K24,K(2,4) 

.(K25,K(2,5) 


),<K34,K<3,4)) 

),(K35,ll(3,5)) 


,<K44,K(4,4> 

,(K45,K(4,5) 


>,(K54,K(5,4>) 

),(K55,K(5,5)> 


,(L21,L<2,1) 

,(L22,L(2,2) 

,(L23,L(2,3) 

,(L24,L(2,4) 

,(L25,L(2,5) 


),(L31,L(3,1)) 

),<L32,L(3,2)) 

),(L33,L(3,3)) 

),{L34,L<3,4)) 

),(L35,L<3,5)) 


,(L41,L(4,1) 

,(L42/L(4,2) 

,(L43,L(4,3) 

,(L44/L(4,4) 

,(L45/L(4,5) 


),(L51,L(5,1)> 

),(L52,L(5,2)) 

),(L53,L(5,3)) 

),(L54,L(5,4)) 

),(L55,L(5,5)) 


C 

CF 


AVE1(K1,K2,K3,K4,K5,K6,K7,K8,K9) - 
( FLOAT(Kl+K2+K3+K4+K5+K6+K7+K8+K9)/9.0 

AVE2(K1,K2,K3,K4,K5,K6,K7) - 
« FLOAT(K1+K2+K3+K4+K5+K6+K7)/7.0 

VARl(Al,Ll,L2,L3,L4,L5,L6,L7,Le,L9) - 
« 7.0*(FLOAT(Ll+L2+L3+L4+L5+L6+L7+L8+L9)-Al*Al*9.0)/9.0 

VAR2(A1,L1,L2,L3,L4,L5,L6,L7) - 
& FLOAT(L1+L2+L3+L4+L5+L6+L7)-A1*A1*7.0 


***** Initialization 


ISXO 

ISYO 

ISXl 

ISYl 

ISX2 

ISY2 

ISX3 

ISY3 


ISX 
ISY 
ISXO-1 
ISYO-1 
ISXO-2 
ISYO-2 
ISXO-3 
ISYO-3 


C 

CF 


***** Smoothing process for the region (3,3>- 
DO 60 IY-3,ISY2 
IYM2 - IY-2 
IYP2 - IY+2 
DO 60 IX-3,ISX2 
IXM2 - IX-2 
IXP2 - IX+2 


■(lSX-2,ISY-2) 


KY - 0 

DO 20 JY-IYM2,IYP2 
KY - KY+1 
KX - 0 

DO 20 JX-IXM2,IXP2 
KX - KX+1 
IPD - IP(JX/JY) 
K(KX,KY) - IPD 
L(KX,KY) - IPD*IPD 
20 CONTINUE 


C 

CF 


*** Calculation of mean values and variances for nine masks 


A1 - AVE1(K22,K32,K42,K23,K33,K43,K24,K34,K44) 

A2 - AVE2(K21,K31,K41,K22,K32,K42,K33) 

A3 - AVE2(K41,K51,K32,K42,K52,K33,K43) 

A4 - AVE2(K42,K52,K33/K43,K53,K44,K54) 

A5 - AVE2(K33,K43,K34/K44,K54,K45,K55) 

A6 - AVE2(K33,K24,K34,K44,K25,K35,K45) 

A7 - AVE2(K23,K33,K14,K24,K34,K15,K25) 

A8 - AVE2(K12,K22,K13,K23,K33,K14,K24) 

A9 - AVE2(K11,K21,K12,K22,K32,K23,K33) 

VI - VAR1(A1,L22,L32,L42,L23,L33,L43,L24;L34,L44) 





V2 - VAR2(A2,L21,L31,M1,L22,L32,L42,L33) 

V3 - VAR2(A3,L41,L51,L32,L42,L52,L33,L43) 

V4 - VAR2(A4,L.42,L52,L33,L43,L53,I/44,L54) 

V5 - VAR2(A5,L33,L43,L34,L44,L54,L45/L55) 

V6 - VAR2(A6,L33,L24,L34,L44,L25,L35,M5) 

V7 - VAR2(A7,L23,L33,U4,L24,L34/L15,L25) 

V8 - VAR2(A8,L12,L22,L13,L23,L33,L14,L24) 

V9 - VAR2(A9,L11,L21,L12,L22,L32,L23,L33) 

CF *** SET THE DIRECTION FOR NINE MASKS 
Q(2)-3.14 
Q(3)-3.925 
Q(4)-4.71 
Q(5)-5.495 
Q(6)-0.0 
Q(7)-0.785 
Q(8)-1.57 
Q(9)-2.355 

CF *** Mean gray value of the mask with minimum variance —> JP(IX/IY) 
RMIN - VI 
MI - 1 
DO 40 1-2,9 

IF (RMIN .LE. V(I)) GO TO 40 
RMIN - V(I) 

MI - I 

40 CONTINUE 

IF(MI.E0.6)THEN 
THET0(IX,IY)-0.0 
GO TO 45 

ELSE IF(MI.EQ.7)THEN 
THET0(IX,IY)-0.785 
GO TO 45 

ELSE IF(MI.E0.8)THEN 
THET0(IX,IY)-1.57 
GO TO 45 

ELSE IF(MI.EQ.9)THEN 
THET0(IX,IY)-2.355 
GO TO 45 

ELSE IF(MI.E0.2)THEN 
THET0(IX,IY)-3.14 
GO TO 45 

ELSE IF(MI.EQ.3)THEN 
THET0(IX,lY)-3.925 
GO TO 45 

ELSE IF(MI.EQ.4)THEN 
THET0(IX,IY)-4.71 
GO TO 45 

ELSE IF(MI.EQ.5)THEN 
THET0(IX,IY)-5.495 
GO TO 45 

ELSE IF(MI.EQ.1>THEN 
GO TO 41 
END IF 

41 CONTINUE 
SMALL-V(2) 

ASMALL-Q(2) 

DO 42 1-3,9 
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IF(SMALL.GT.V(1))THEN 
SMALJ>V( I) 


ASMALL-O(I) 

ENDIF 
42 CONTINUE 

THETO(IX,IY)-ASMALL 
45 CONTINUE 

JP(IX,IY) - A(MI)+0.5 


60 CONTINUE 
C 

CF ***** Operation for the four corners 
C 

C- 

C 

DO 80 IY-1;4 
DO 80 IX-1,4 
IPD - IP(IX,1Y) 

K(IX,IY) - IPD 
L(IX,IY) - IPD*IPD 
80 CONTINUE 
C 

CF***** JP(1,1) ***** 

JP(1,1) - AVE2(Kll,K21,K12,K22,K32,K23,K33)+0.5 
C 

CF***** JP(2,1) ***** 

A1 - AVE2(K21,K31,K22,K32,K42,K33,K43) 

A2 - AVE2(K21,K12,K22,K32,K13,K23,K33) 

VI - VAR2(A1,L21,L31,L22,L32,L42,L33,L43) 

V2 - VAR2(A2,L21,L12,L22,L32,L13,L23,L33) 

C 

RMIN - A1 

IF (VI .GT. V2) RMIN - A2 
JP(2,1) - RMIN+0.5 
C 

CF***** JP(1,2) ***** 

A1 - AVE2(K21,K31,K12,K22,K32,K23,K33) 

A2 - AVE2(K12,K22,K13,K23,K33,K24,K34) 

VI - VAR2(A1,L21,L31,L12,L22,L32,L23,L33) 

V2 - VAR2(A2,L12,L22,L13,L23,L33,L24,L34) 

C 

RMIN - A1 

IF (VI .GT. V2) RMIN - A2 
JP(1,2) - RMIN+0,5 
C 

CF***** JP(2,2) ***** 

A1 - AVE1(K11,K21,K31,K12,K22,K32,K13,K23,K33) 

A2 - AVE2(K31,K41,K22,K32,K42,K33,K43> 

A3 - AVE2(K22,K32,K23,K33,K43,K34,K44) 

A4 - AVE2(K22,K13,K23,K33,K14,K24,K34) 

VI - VAR1(A1,L11,L21,L31,L12,L22,L32,L13,L23,L33) 
V2 - VAR2(A2,L31,L41,L22,L32,L42,L33,L43) 

V3 - VAR2(A3,L22,L32,L23,L33,L43,L34,L44) 

V4 - VAR2(A4,L22,L13,L23,L33,L14,L24,L34) 

C 

RMIN - VI 
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non 


MI - 1 

DO 100 1-2,4 

IF (RMIN .LE. V(I)) GO TO 100 
RMIN - V<I) 

MI - I 

100 CONTINUE 

JP{2,2) - A(MI)+0.5 


DO 120 n-1/4 
KX - 0 

DO 120 IX-ISX3,ISX0 

KX - KX+1 

IPD - 1P(IX,IY) 

K(KX,IY) - IPD 
L(KX,IY) - IPD*IPD 
120 CONTINUE 
C 

CF***** JP(ISX-1,1) ***** 

A1 - AVE2<K31,K22,K32,K42,K23,K33#K43) 

A2 - AVE2(K21,K31,K12,K22,K32,Kp»X23) 

VI - VAR2(A1,L31,L22,L32,M2,L23,L33,W3) 

V2 - VAR2(A2,Ii21,L31,L12,L22,L32,L13,L23) 

C 

RMIN - A1 

IF (VI .GT. V2) RMIN - A2 
JP{ISX1,1) - RMIN+0.5 

c 

rr***** JP(ISX”1,2) ***** 

A1 - AVE1(K21,K31,K41,K22,R32,K42,K23,K33,K43) 

A2 - AVE2(K32,K23,K33,K43,K24,K34,K44) 

A3 - AVE2(K22,K32,K13,K23,K33,K14,K24) 

A4 - AVE2(K11,K21,K12,K22,K32/K13/K23) 

VI - VAR1(A1,L21,L31,L41,L22,L32,L42,L23,L33,L43) 
V2 - VAR2(A2,L32,L23,L33,L43,L24,L34,M4) 

V3 - VAR2(A3,L22,L32,L13,L23,L33,L14,L24) 

V4 - VAR2{A4,L11,L21,L12,L22,L32,L13,L23) 

C 

RMIN - VI 
Ml - I 

DO 140 1-2,4 

IF (RMIN .LE. V(I)) GO TO 140 
RMIN - V(I) 

MI - I 

140 CONTINUE 
C 

JP(ISX1,2) - A(Ml)+0.5 

-’’aTC2;K32,M2!k23.1!33,K43,K24,M4) 

A2 - AVE2(K21,K31,K22,K32,K42,K23,Kp) 

VI - VAR2(A1,L32,L42,L23,L33,L43,L24,L34) 

V2 - VAR2(A2,L21,L31,L22,L32,L42,L23,L33) 



uuo 


RMIN - A1 

IF (VI .GT. V2) RMIN - A2 
JP(ISX0,2> - RMIN+0.5 


KY - 0 

DO 160 IY-ISY3,ISY0 
KY - KY+1 
DO 160 IX-1,4 
IPD - 1P(IX,1Y) 

K(IX,KY) - IPD 
L(IX,KY) - IPDMPD 
160 CONTINUE 
C 

CF***** JP(1,ISY-1) ***** 

A1 - AVE2(K21,K31,K12,K22,K32,K13,K23) 

A2 - AVE2(K22,K32,K13,K23,K33,K24,K34). 

VI - VAR2(A1,L21/L31,L12,L22,L32,L13,L23) 

V2 - VAR2(A2,L22,L32,L13,L23,L33,L24,L34) 

C 

RMIN - A1 

IF (VI .GT. V2) RMIN - A2 
JP(1,ISY1) - RMIN+0.5 
C 

CT***** JP(2,ISY-1) ***** 

A1 - AVE1(K12,K22,K32,K13,K23,K33,K14,K24,K34) 

A2 - AVE2(K11,K21,K31,K12,K22,K32,K23) 

A3 - AVE2(K31,K41,K22,K32,K42,K23,K33) 

A4 - AVE2(K32,K42,K23,K33,K43,K34,K44) 

VI - VAR1(A1,L12,L22,L32,L13,L23,L33/L14,L24,L34). 
V2 - VAR2(A2/L11,L21,L31,L12,L22,L32,L23> 

V3 - VAR2(A3,L31,L41,L22,L32,L42,L23,L33) 

V4 - VAR2(A4,L32,L42,L23,L33,L43,L34,L44) 

C 

RMIN - VI 
MI - 1 
DO 180 1-2,4 

IF (RMIN .LE. V(I)) GO TO 180 
RMIN - V(I) 

MI - I 

180 CONTINUE 
C 

JP(2,ISY1) - A(MI)+0.5 
C 

CF***** JP(1,ISY) ***** 

JP(1,ISY) - AVE2(K22,K32,K13,K23,K33,K14,K24)+0.5 
C 

CF***** ■JP(2,ISY) ***** 

A1 - AVE2(K12,K22,K32,K13,K23,K33,K24) 

A2 - AVE2(K32,K42,K23,K33,K43,K24,K34) 

VI - VAR2(A1,L12,L22,L32,L13,L23,L33,L24) 

V2 - VAR2(A2,L32,L42,L23,L33,L43,L24,L34) 

C 

RMIN - A1 

IF (VI .GT. V2) RMIN - A2 
JP(2,ISY0) - RMIN+0.5 
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KY - 0 

DO 200 IY-1SY3,ISY0 
KY - KY+1 
KX - 0 

DO 200 1X-ISX3,ISX0 
' KX - KX+1 

IPD - IP(IX,IY) 

K(KX,KY) - IPD 
L(KX,KY) - IPD*IPD 
200 CONTINUE 
C 

CF***** JP(ISX-1,1SY-1) ***** 

A1 - AVE1(K22,K32,K42,K23,K33,K43,K24,K34,K44) 

A2 - AVE2(K21,K31,K41,K22,K32,K42,K33) 

A3 - AVE2(K12,K22,K13,K23,K33,K14,K24) 

A4 - AVE2(K11,K21,K12,K22,K32,K23,K33) 

VI - VAR1(A1,L22,L32,L42,L23,L33,L43,L24,L34,L44) 

V2 - VAR2(A2,L21,L31,L41,L22,L32,L42,L33) 

V3 - VAR2(A3,L12,L22,L13,L23,L33,L14,L24) 

V4 - VAR2(A4,L11,L21,L12,L22,L32,L23,L33) 

C 

RMIN - VI 
MI - I 
DO 220 1-2,4 

IF (RMIN .LE. V(l)) GO TO 220 
RMIN - V(I) 

MI - 1 

220 CONTINUE 
C 

JP(ISX1,ISY1) - A(MI)+0.5 
C 

CF***** JP(ISX,ISY-1) ***** 

A1 - AVE2(K22,K32,K23,K33,K43,K24,K34) 

A2 - AVE2(K21,K31,K22,K32,K42,K33,K43) 

VI - VAR2(A1,L22,L32,L23,L33,L43,L24,L34) 

V2 - VAR2(A2,L21,L31,L22,L32,L42,L33,L43) 

C 

RMIN - A1 

IF (VI .GT. V2) RMIN - A2 
JP(ISX0,ISY1) - RMIN+0.5 
C 

CF***** JP(ISX-1,ISY) ***** 

A1 - AVE2(K22,K32,K42,K23,K33,K43,K34) 

A2 - AVE2(K12,K22,K13,K23,K33,K24,K34) 

VI - VAR2(A1,L22,L32,L42,L23,L33,L43,L34) 

V2 - VAR2(A2,L12,L22,L13,L23,L33,L24,L34) 

C 

RMIN - A1 

IF (VI .GT. V2) RMIN - A2 
JP(ISX1,ISY0) - RMIN+0.5 
C 

CF***** JP(1SX,ISY) ***** 

JP(ISX0,ISY0) - AVE2(K22,K32,K23,K33,K43,K34,K44)+0.5 
C 

CF ***** Marginal operation 
C 

C- 

C 
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DO 360 1X-3,ISX2 
IXM2 - IX-2 
IXP2 - IX+2 
C 

Q - 

C 

DO 240 JY-1,4 
KX - 0 

DO 240 JX-IXM2,IXP2 

KX - KX+1 

IPD - 1P(JX,JY) 

K(KX,JY) - IPD 
L(KX,JY) - IPD*IPD 
240 CONTINUE 
C 

CF***** JP(1X,1) ***** 

A1 - AVE2(K31,K41,K32,K42,K52,K43,K53) 

A2 - AVE2(K31,K22,K32,K42,K23,K33,K43) 

A3 - AVE2(K21,K31,K12,K22,K32,K13,K23) . 

VI - VAR2(A1,L31,I/41,L32,L42,L52,L43,L53) 

V2 - VAR2<A2,L31,L22,L32,L42,L23,L33,L43) 

V3 - VAR2(A3,L21,L31,L12,L22,L32,L13,L23) 

C 

RMIN - VI 
MI - 1 
DO 260 1-2,3 

IF (RMIN .LE. V(I)) GO TO 260 
RMIN - V(I) 

MI - I 

260 CONTINUE 
C 

JP(IX,1) - A(MI)+0.5 
C 

CF***** JP(IX,2) ***** 

A1 - AVE1(K21,K31,K41,K22,K32,K42,K23,K33,K43) 

A2 - AVE2(K41,K51,K32,K42,K52,K43,K53) 

A3 - AVE2(K32,K42,K33,K43,K53,K44,K54) 

A4 - AVE2(K32,K23,K33,K43,K24,K34,K44) 

A5 - AVE2(K22,K32,K13,K23,K33,K14,F24> 

A6 - AVE2(K11,K21,K12,K22,K32,K13,K23) 

VI - VAR1(A1,L21,L31,L41,L22,L32,L42,L23,L33,L43) 
V2 - VAR2(A2,L41,L51,L32,L42,L52,L43,L53) 

V3 - VAR2(A3,L32,L42,L33,L43,L53,L44,L54) 

V4 - VAR2(A4,L32,L23,L33,L43,L24,L34,L44) 

V5 - VAR2(A5,L22,L32,L13,L23,L33,L14,L24) 

V6 - VAR2(A6,L11,L21,L12,L22,L32,L13,L23) 

C 

RMIN - VI 
MI - 1 
DO 280 1-2,6 

IF (RMIN .LE. V(I)) GO TO 280 
RMIN - V(I) 

MI - I 

280 CONTINUE 
C 

JP(IX,2) - A(MI)+0.5 
C 

C- 
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C 

KY - 0 

DO 300 JY-ISY3,ISY0 
KY - KY+1 
KX - 0 

DO 300 JX-IXM2,IXP2 ‘ 

KX - KX+1 
IPD - IP(JX,JY) 

K(KX,KY) - IPD 
L(KX,KY) - IPD*IPD 
300 CONTINUE 
C 

CT***** JP(IX,ISY-1) ***** 

M - AVE1(K22,K32,K42,K23,K33,K43,K24,K34,K44) 

A2 - AVE2(K42,KS2,K33«K43,K53,K44,KS4) 

A3 - AVE2(K41,K51,K32,K42,K52,K33,K43) 

A4 - AVE2(K21,K31,K41,K22,K32,K42,K33) 

A5 - AVE2(K11,K21,K12/K22,K32,K23,K33) 

A6 - AVE2(K12,K22,K13/K23,K33,K14,K24) 

VI - VARl(Al,L22,li32/L42,L23,L33,l/43,L24,L34,L44> 
V2 - VAR2(A2»L42,L52,L33,L43,L53«L44,L54) 

V3 - VAR2(A3,L41,L51,L32,L42,L52,L33,L43) 

V4 - VAR2(A4,L21,L31,L41,L22,L32,L42,L33) 

V5 - VAR2(A5,L11,L21,L12,L22,L32,L23,L33) 

V6 - VAR2(A6,L12,L22,L13,L23,L33,L14,L24) 

C 

RMIN - VI 
MI - 1 
DO 320 1-2,6 

IF (RMIN .LE. V(I)) GO TO 320 
RMIN - V(I) 

MI - I 

320 CONTINUE 
C 

JP(IX,ISY1) - A(MI)+0.5 
C 

QY***** JP(IX,ISY) ***** 

A1 - AVE2(K42,K52,K33,K43,K53,K34,K44) 

A2 - AVE2(K22,K32,K42,K23,K33,K43,K34) 

A3 - AVE2(K12,K22,K13,K23,K33,K24,K34) 

VI - VAR2(A1,L42,L52,L33,L43,L53,L34,L44) 

V2 - VAR2(A2,L22,L32,L42,L23,L33,L43,L34) 

V3 - VAR2{A3,L12,L22,L13,L23,L33,L24,L34> 

C 

RMIN - VI 
MI - 1 
DO 340 1-2,3 

IF (RMIN .LE. V(I)) GO TO 340 
RMIN - V(I) 

MI - I 

340 CONTINUE 
C 

JP(IX,ISY0) - A(MI)+0.5 
C 

360 CONTINUE 
C 
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DO 500 1Y-3,ISY2 
IYM2 - lY-2 
IYP2 - IY+2 
C 

Q - 

C 

KY - 0 

DO 380 JY-IYM2,IYP2 
KY - KY+1 
DO 380 JX-1,4 
IPD - IP(JX,JY) 

K(JX,Ky) - IPD 
L(JX,KY) - IPD*IPD 
380 CONTINUE 
C 

CF****« JP(1,IY) ***** 

A1 - AVE2(K21,K31,K12,K22,K32,K13,K23) 

A2 - AVE2(K22,K32,K13,K23,K33,K24,K34> 

A3 - AVE2(K13,K23,K14,K24,K34,K25,K35) 

VI - VAR2(A1,L21,L31,L12,L22,L32,L13,L23) 

V2 - VAR2(A2,L22,L32,L13,L23/L33,L24,L34) 

V3 - VAR2(A3,L13,L23,L14,L24,L34,L25,L35) 

C 

RMIN - VI 
MI - 1 
DO 400 1-2,3 

IF (RMIN .Lt. V(I)) GO TO 400 
RMIN - V(I) 

MI - I 

400 CONTINUE 
C 

JP(1,IY) - A(MI)+0.5 
C 

CF***** JP(2,IY) ***** 

A1 - AVE1(K12,K22,K32,K13,K23,K33,K14,K24,K34) 

A2 - AVE2(K11,K21,K31,K12,K22,K32,K23) 

A3 - AVE2(K31,K41,K22,K32,K42,K23,K33) 

A4 - AVE2(K32,K42,K23,K33,K43,K34,K44) 

A5 - AVE2(K23,K33,K24,K34,K44,K35,M5) 

A6 - AVE2(K23,K14,K24,K34,K15,K25,K35) 

VI - VAR1(A1,L12,L22,L32,L13,L23,L33,L14,L24,L34) 
V2 - VAR2(A2,L11,L21,L31,L12,L22,L32,L23) 

V3 - VAR2(A3,L31,L41,L22,L32,L42,L23,L33) 

V4 - VAR2(A4,L32,L42,L23,L33,L43,L34,L44) 

V5 - VAR2(A5,L23,L33,L24,L34,L44,L35,L45) 

V6 - VAR2(A6,L23,L14,L24,L34,L15,L25,L35) 

C 

RMIN - VI 
MI - 1 
DO 420 1-2,6 

IF (RMIN .LE. V(I)) GO TO 420 
RMIN - V(I) 

MI - I 

420 CONTINUE 
C 

JP(2,iy) - A(MI)+0.5 
C 

C- 
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c 

KY - 0 

DO 440 JY-IYM2,IYP2 
KY - KY+1 
KX - 0 

DO 440 JX-ISX3,ISX0 

KX - KX+1 

IPD - IP{JX,JY) 

K(KX,KY) - IPD 
L(KX,KY) - IPD*1PD 
440 CONTINUE 
C 

CF***** JP{ISX-1,IY) ***** 

A1 - AVE1(K22,K32,K42,K23,K33,K43,K24,K34,K44^ 

A2 - AVE2(K21«K31,K41,R22,K32,K42,K33) 

A3 - AVE2(K11,K21,K12,K22,K32,K23,K33) 

A4 - AVE2(K12,K22/K13,K23,X33/K14,X24) 

A5 - AVE2(K23,K33,K14,K24,K34,K15,K25) 

A6 - AVE2(K33,K24,K34,K44,K25,K35,K45) 

VI - VAR1{A1,L22,L32,L42,L23,L33,L43,L24,L34,L44) 
V2 - VAR2(A2,L21,L31,L41,L22,L32,L42,L33) 

V3 - VAR2(A3,L11,L21,L12,L22,L32,L23,L33) 

V4 - VAR2(A4,L12,L22,L13,L23,L33,L14,L24) 

V5 - VAR2(A5,L23,L33,L14,L24,L34,L15,L25) 

V6 - VAR2(A6,L33,L24,L34,L44,L25,L35,L45) 

C 

RMIN - VI 
MI - 1 
DO 460 1-2,6 

IF (RMIN .LE. V(I)) GO TO 460 
RMIN - V(I) 

MI - I 

460 CONTINUE 
C 

JP(ISX1,IY) - A(MI)+0.5 
C 

CF***** JP(ISX,1Y) ***** 

A1 - AVE2(K21,K31,K22,K32,K42,K33,K43) 

A2 - AVE2(K22,K32,K23,K33,K43,K24,K34) 

A3 - AVE2(K33,K43,K24,K34,K44,K25,K35) 

VI - VAR2(Al,L2l,L31,L22,L32,L42,L33,L43) 

V2 - VAR2(A2,L22,L32,L23,L33,L43,L24,L34) 

V3 - VAR2(A3,L33,L43,L24,L34,L44,L25,L35) 

C 

RMIN - VI 
MI - 1 
DO 480 1-2,3 

IF (RMIN .LE. V(I)) GO TO 480 
RMIN - V(I) 

MI - I 

480 CONTINUE 
C 

JP(ISX0,IY) - A(MI)+0.5 
C 

500 CONTINUE 

DO 600 1-1,256 
DO 600 J-1,256 
IP(I,J)-JP(I,J) 

600 CONTINUE 
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DO 601 IX-1,2 
DO 601 IY-1,256 
THET0(IX,IY)-0.0 

601 CONTINUE 

DO 602 IX-1,256 
DO 602 IY-1,2 

THETO(IX,IY)-0.0 

602 CONTINUE 

DO 603 IX-255,256 
DO 603 IY-1,256 
THET0(IX,IY)-0.0 

603 CONTINUE 

DO 604 IX-1,256 
DO 604 IY-255,256 
THETO(IX,IY)-0.0 

604 CONTINUE 

DO 605 IX-1,256 
DO 605 IY-1,256 

THETl(IX,lY)-THETO (IX,lY) 

605 CONTINUE 

DO 606 IX-1,256 
DO 606 IY-1,256 

IF{THET1(IX,lY) .GE. 3.14)THEN 

THETl(IX,lY)-THETl(IX,IY)-3.14 
ELSE 

THETl(IX,lY)-THETl(IX,lY) 

ENDIF 

606 CONTINUE 
C 

CF ***** RETURN 
RETURN 
END 

SUBROUTINE RXEP(IP,THETl,PRBEl,ISX,ISY,TEMP) 

C 

C Copyright (c) 1983 by AIST MITI(JAPAN) 

C 

CP Computes initial parameter and initial probability for edge 
CP reinforcement by relaxation method (Lower-level routine 

CP for subroutine RXEG). 

C 

CS CALL RXEP(IP,THETl,PRBEl,ISX,ISY) 

C 

CK RELAXATION, EDGE 

C 

CA IP(ISX,ISY) ! Input image 

CA THETl(ISX,ISY): Initial value of gradient direction for 

CA each pixel 

CA PRBE1(ISX,ISY)s Initial value of edge probability for 

CA each pixel 

C 

CN Reference 

CN @1! B. J. Schachter, A. Lev, S. W. Zucker, and A. Rosenfeld, "An 
CN application of relaxation methods to edge reinforcement," 

CN IEEE Trans. Syst., Man, Cybern., vol. SMC-7, pp. 813-816, 

CN Nov. 1977. 

C 

CD JULY 1979. PROGRAMMED BY K.SAKAUE 

C 

CM A small real number -1E30 is set in variable AMX as an initial 
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CM value. When using a computer with different word lengths/ care 

CM should be taken. 

C 

DIMENSION 1P(ISX,ISY)/TEMP(ISX,ISY) 

DIMENSION THET1(ISX,ISY) 

DIMENSION PRBE1(ISX/ISY) 

DIMENSION NE(3/3) 

C 

CF INITIAL EDGE VALUES 

CF BY PREWITT OPERATOR (DIFFERENTIAL TYPE) 

C 

AMX—1E30 
DO 10 IY-1,ISY 
DO 10 IX-1,ISX 
DO 11 J-1,3 
DO 11 1-1,3 
Jl-IY+J-2 
ll-IX+I-2 

IF(J1.LT.1.0R.J1.GT.ISY)G0 TO 12 
IF(I1.LT.1.0R.I1.GT.ISX)G0 TO 12 
NE(I,J)-1P(11/J1) 

11 CONTINUE 

DX-FLOAT(NE(3/l)+NE(3,2)+NE(3,3)-NE(l,l)-NE(l/2)-NE(l/3)) 
DY-FLOAT(NE(1/3)+NE(2,3)+NE(3,3)-NE( 1,1)-NE(2/1)-NE(3,1)) 
AMG-SQRT(DX* * 2+DY * * 2) 

IF(AMG.GT.AMX)AMX-AMG 
PRBE1(1X/IY)-AMG 

IF(DX.LE.0E0.AND.DY.LE.OEO)DX-0.0001 


GO TO 10 
12 CONTINUE 

PRBE1(IX/IY)-0.001 


C 

C 

C 

CS 

C 

CP 

CP 

CP 

C 

CK 

C 

CA 

CA 

CA 

CA 


10 CONTINUE 

DO 20 IY-1,ISY 
DO 20 IX-1,ISX 

PRBEl(IX,IY)-PRBE1(IX,IY)/AMX 
PRBE1(IX,IY)-AMAX1(.01,AMIN1(.9,PRBE1(IX,lY))) 

TIMP(IX,lY)-PRBEl(IX,lY) 

20 CONTINUE 
RETURN 
END 

SUBROUTINE RXEI (THETl, THET2, PRBEl, PRBE2, ISX, ISY, ISNX, ISNY, C, W) 


Copyright (c) 1983 by AIST MITI(JAPAN) 


CALL RXEI(THETl,THET2,PRBEl,PRBE2,ISX,ISY,ISNX,ISNY,C,W) 

Updates edge parameter and edge probed>ility for edge 
reinforcement by relaxation method. 

(Lower-level routine for subroutine RXEG) 


RELAXATION, EDGE 

THET1(ISX,ISY): 
THET2(ISX,ISY)! 
PRBE1(ISX,ISY): 
PRBE2(ISX,ISY)! 


Gradient direction 
Gradient direction 
Edge probabilities 
Edge probabilities 


(before updating) 
(after updating) 
(before updating) 
(after updating) 
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C 


c 

c 

c 


c 

CF 

C 


C 

CF 

CF 

c 


ISNX/ISNY ! Size of neighborhood being considered 

(limited to odd numl'er three or higher) 
C(4) : Parameter Cl through C4 (See the 

reference.) 

.866,.124,.005,.005 (standard) 

.706, .176, .059 (noise cleaning) 

W : Parameter W (See the reference.) 

3.0 for standard use. 

DIMENSION THET1(ISX,ISY) 

DIMENSION THET2(ISX,-TSY) 

DIMENSION PRBE1(1SX,ISY) 

DIMENSION PRBE2(ISX,ISY) 

DIMENSION C(4) 

DATA HPAI/1.570796327/ 


ISFTX-ISNX/2+1 
ISFTY-ISNY/2+1 
DO 10 IY-1,ISY 
DO 10 IX-1,ISX 

ALPHA-THET1(IX,IY) 

PXY-PRBE1(IX,lY) 

Ql-0.0 

Q2-0.0 

DHX-W* PXY*C0S(ALPHA) 

DHY-W*PXY*SIN(ALPHA) 

DO 11 J-1,ISNY 
DO 11 I-1,ISNX 
lU-IX+l-ISFTX 
IV-IY+J-ISFTY 

IF(rV.EQ.lY.AND.IU.EO.IX)GO TO 11 

IF(IV.LT.1.0R.IV.GT.ISY)G0 TO 11 

IF(IU.LT.1.0R.IU.GT.ISX)G0 TO 11 

IDY-IV-IY 

IDX-IU-IX 

JI>-IABS(IDY) 

ID-IABS(IDX) 

LD-MAX0(ID,JD) 

LD ! CHESSBOARD DISTANCE FROM (IX,IY) TO (IU,R') 

ALD-FLOAT(LD) 

D2-1.0/2.0**ALD 

BETA-THET1(IU,IV) 

DY-FLOAT(IDY) 

DX-FLOAT(IDX) 

IF(IDX.EQ.0.AND.IDY.EQ.0)DX-.00001 
GAMMA«ATAN2(DY,DX) 

PUV-PRBE1(IU,IV) 

COMPATIBILITY COEFFICIENTS 
EDGE/EDGE INTERACTION 

IF(ALPHA.EQ.BETA.AND.BETA.EQ.GAMMA)THEN 
REE-1.0 


ELSE IF(ALPHA. EQ. BETA. AND. GAMMA, NE. BETA) THEN 
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REE-0.5 
ELSE 

REE-COS (ALPHA-GAMMA) *COS (BETA-GAMMA) * D2 
ENDir 

c 

CF INTERACTION WITH NONEDGES 

C 

AG2-2.0*(ALPHA-GAMMA) 

AG2—COS(AG2)*D2 
REN-AMIN1(0.0,AG2) 

BG2-2.0*(BETA-GAMMA) 

RNE-(1.O-COS(BG2))*D2*0.5 
RNN-D2 
C 

CF COMBINED REINFORCEMENT PROCESS 

C 

Q1-Q1+C(1)*PUV*REE+C(2)*(1.0-PUV)*REN 
02-Q2+C(3)*PUV*RNE+C(4)*(1.0-PUV)*RNN 
PREE-PUV*REE 
DHX-DHX+PREE*COS(BETA) 

DHY-DHY+PREE* SIN(BETA) 

11 CONTINUE 

00-ABS(Q1)+ABS(02) 

Ql-Ql/QQ 

Q2-Q2/QO 

P1D-PXY*(1.0+Q1) 

P2D-(1.0-PXY)*(1.0+Q2) 

PNEW-PlD/(P1D+P2D) 

IF(DHX.LE.OEO.AND.DHY.LE.OEO)DHY-0.00001 
THNEW-ATAN2(DHY,DHX) 

PRBE2(IX,IY)-PNEW 
THET2(IX,1Y)-THNEW 
10 CONTINUE 
RETURN 
END 
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